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We demonstrate how the spin-1/2 XYZ quantum Heisenberg model can be realized with bosonic 
atoms loaded in the p-band of an optical lattice in the Mott regime. The sign and relative strength 
of the couplings characterizing the model are demonstrated to be experimentally tuneable. We 
discuss the phase diagram in the one dimensional case, and show that finite size effects relevant for 
trapped systems lead to devil's staircase structure. Finally, we discuss experimental issues related 
to preparation, manipulation and detection. 
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Introduction. — Powerful tools developed recently to 
unravel the physics of many-body quantum systems offer 
an exciting platform for understanding quantum mag- 
netism. It is now possible to engineer different systems 
in the lab that mimic the physics of theoretically chal- 
lenging spin models [l|, thereby performing "quantum 
simulations" as suggested long ago [2|. Along these lines, 
systems of trapped ions are promising candidates and 
have already been employed to simulate both small [3[ 
and large [j| numbers of spins. In these setups, however, 
sustaining control over the parameters becomes very dif- 
ficult as the system size increases. Furthermore, due to 
confining trapping potentials, realizations are limited to 
chains with up to 25 spins. It is also very difficult to 
construct paradigmatic spin models with short range in- 
teractions using systems of trapped ions. 

Such obstacles are less severe for realizations of spin 
models with cold atoms in optical lattices [l|. A bosonic 
system in a deep tilted lattice has recently been used to 
simulate the phase transition in the ID Ising model [5[- 
Fermionic atoms were employed to study dynamical 
properties of quantum magnetism beyond spin-1/2 sys- 
tems |6| , in which case the spin degrees of freedom were 
encoded in the internal atomic states. This idea, first in- 
troduced in Ref . [7| , has also been applied to other config- 
urations, and simulation of different types of spin models 
have been proposed |8|. However, due to the character of 
the atomic s-wave scattering among the different Zeeman 
levels, such mappings usually yield effective spin models 
supporting continuous symmetries like the XXZ model, 
which is invariant under U{1) transformation. But as the 
main goal of a quantum simulator is to realize systems 
that cannot be tackled via analytical and/or numerical 
approaches, it is important to explore alternative scenar- 
ios that yield low symmetry spin models with anisotropic 
couplings and external fields. 

In this paper we propose such a scenario by demon- 
strating that bosonic atoms in the first excited band 
(p-band) of a two-dimensional (2D) optical lattice can 



realize the spin-1/2 XYZ quantum Heisenberg model 
in an external field. Systems of cold atoms in excited 
bands feature an additional orbital degree of freedom [9| 
that gives rise to novel physical prop erties whose theo- 
retical [10| and experimental [llj, |l2| aspects have been 
recently investigated. Their dynamics include "orbital 
changing" interactions, where two atoms in one of the 
orbital states scatter into two atoms that occupy orbital 
states of a different type. This is the key mechanism 
leading to the anisotropy of the effective spin model: 
These processes reduce the continuous t/(l) symmetry 
into a set of discrete Z2 symmetries characteristic of the 
XYZ model. We furthermore show how the strength 
and sign of the couplings can be experimentally con- 
trolled. This means that one can realize a whole class 
of anisotropic XYZ models with ferromagnetic and/or 
anti-ferromagnetic correlations. 



To illustrate the rich physics which can be explored 
with this system and to examine finite size effects rele- 
vant for the trapped systems, we discuss the phase dia- 
gram for the ID XYZ chain. This case exhibits ferromag- 
netic as well as anti-ferromagnetic phases, a magnetized 
field-dependent state or polarized phase, a spin-fiop, and 
a fioating phase jl3[. We provide an exact diagonaliza- 
tion of this model for a finite system relevant for the 
trapped case. This reveals that finite size effects result 
in the appearance of a devil's staircase manifested in the 
form of spin density waves. Finally, we discuss how to 
experimentally probe and manipulate the spin degrees of 
freedom. 



p- orbital Base system. — We consider bosonic atoms 
of mass m in an optical lattice of the form V{r) = 
Vx sin^ {k^x) + VySin^{kyy) + VzSm'^{kzz) with Va and 
ka the amplitude and wave number in the direction a. 
Assuming that all atoms arc in the first excited p orbital. 



the tight-binding Hamiltonian is 
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Here a] ^ creates a bosonic particle in the orbital a = 

Px,Py,Pz at site i, hi^a = CLJ^ai^a, and the sum is over 
nearest neighbors i, j. The tunneling matrix elements are 
given by tf^ = -Jdrwf{r)* [-h^V'^ /2m + V{r)] wf{r) 
where wf (r) is the Wannier function of orbital a at site 
i. Note that tf is anisotropic. For instance, a boson in 
the Pa;-orbital has a much larger tunneling rate in the x- 
direction than in the j/-direction. The coupling constants 
are given by Uaa' = Uo J dr \wf {r)\'^ \wf {r)\'^ , with Uq 
(> 0) the onsite interaction strength determined by the 
scattering length. The last term in ([T]) is the orbital 
changing term describing the flipping of a pair of atoms 
from the state a' to the state a. 

Effective spin Hamiltonian. — We are interested in the 
physics of the Mott insulator phase with unit filling in the 
strongly repulsive limit \tfj'^ <C Uaa'- Projecting onto 
the Mott space of singly occupied sites with the opera- 
tor P, the Schrodinger equation becomes ^Mott-PlV^) = 
EP\iP) with i^Mott == -PH{Hq - E)-^HP. Here Q = 
1 - P and iJg = QHQ [3 [H. Since we will find that 
E - t^/U, we can take {Hq - E)''^ = H'^ . 

For 2D, the space of doubly occupied states of a given 
site j is three-dimensional and spanned by \pxPx) = 
2-'/'sl«]j0), \pypy) = 2-i/2a]^a]j0), and \pxPy) ^^ 
dj^aj |0). In this space, it is straightforward to find Hq 
from ([T]), and subsequent inversion yields 



Uyy/U^ -Uxy/U^ 






(2) 



with U^ = UxxUyy — U^y. Using ^ we can now calculate 
all possible matrix elements of Hmou hi the Mott space. 
After some algebra we find 
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map ([3]) onto a spin-1/2 XYZ model in an external field 
HxYZ = Y. J., [(1 + 7)Sr^; + (1 - 7)^?5J 
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Here, {i,j) means summing over each nearest neighbor 
pair i,j only once. The couplings are given by Jij — 
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where Ea is the onsite energy of the orbital a. Note that 
t^jt^ji < for p-band bosons [9| and therefore Jij > 0. 
Also since I7I < 1, the interactions between the y- 
component of neighboring spins favor anti-ferromagnetic 
order. As we shall see below, the sign of the couplings 
for interactions between the spin components a;- and be- 
tween Z-, as well as the relative strength of all couplings 
can be experimentally controlled. 

Equation ([4]) is a main result of this paper. It demon- 
strates how p-orbital bosons in a 2D optical lattice can re- 
alize the XYZ quantum spin-1/2 Heiscnberg model. The 
fact that 7 7^ is a consequence of the orbital changing 
term in ([T]), and it implies that the continuous U{1) sym- 
metry for S^ and S^ is broken down to a set of Z2 sym- 
metries. The Z2 symmetries reflect the 'parity' conser- 
vation in the original bosonic picture which classifies the 
many-body states according to total even or odd number 
of atoms in the px- and p^-orbitals. We emphasize that 
this derivation makes no assumptions regarding the ge- 
ometry of the 2D lattice - i.e. it can be square, hexagonal 
etc. 

To illustrate the rich physics of the XYZ model and to 
examine finite size effects, we now focus on the case of a 
ID lattice where we can obtain a reliable numerical solu- 
tion of (U, and where quantum fluctuations are strong. 
Note that by increasing both the lattice amplitude and 
spacing in the y-direction keeping Vyky c^Vxk"^, one can 
exponentially suppress tunneling in the y-direction to ob- 
tain a ID model, while the px and py orbitals arc still 
quasi-degenerate [l6|. In the ID setting, we will drop 
the "direction" subscript ij on the coupling constants. 

For ID, the importance of the orbital changing term 
can be further illuminated by employing the Jordan- 
Wigner transformation S^ = e*'^^j=i '^^'^^ Ci for fermionic 
operators Ci to the spin chain. The result is the fermionic 
Hamiltonian 



(3) Hk/J - Y. [(4c«+i 
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where x = y, and y = x. By further employ- 
ing the Schwinger angular momentum representation, 

S! = Wl^axi - al^dyi), S+ = Sf + iSf = al^ay, and 
S~ ~ Sf — iSf = a/ -cixi, together with the constraint 
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'^IriAxi+o.l^i&yi — Ij wc Can (ignoring irrelevant constants) 



jicicn + l/2)(cl+iC„+i - 1/2) + j{cicn - 1/2) 



(5) 



We see that 7 7^ leads to a pairing term that typically 
opens a gap in the energy spectrum. Incidentally the 
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FIG. 1. (Color online) (a) Schematic phase diagram of the 
XYZ chain, (b) Finite size 'phase diagram' obtained by exact 
diagonalization of 18 spins. The finite size 'phase diagram' 
comprises an incomplete devil's staircase of SDW between 
the PP and AFM phases. As the system size is increased, the 
devil's staircase approaches the SF and the FP of (a). The 
anisotropy parameter is 7 = 0.2 in (b). 



limit of A — !■ in Eq. ([5]) is a realization of the Kitaev 
chain (T^j. 

XYZ phase diagram. — We focus here on the most nat- 
ural case for the implementation with p-orbital bosons, 
namely when the x— and y— components of spin favor 
anti- ferromagnetic ordering (I7I < 1). The schematic 
phase diagram is illustrated in Fig. [T] (a). At zero field, 
the XYZ model is integrable 18|. For large positive val- 
ues of A/J the system is anti-ferromagnetic (AFM) in 
the z-direction. Small values of A/J are characterized by 
Neel ordering in the y— direction and the system is in the 
so-called spin- flop phase (SF). The /i = line for large 
negative values of A/J is characterized by a ferromag- 
netic phase (FM) in the z-direction, and for all the cases, 
the limit of large external field displays a magnetized 
phase (PP), where the spins align along the orientation 
of the field in the z-direction. These three phases also 
characterize the phase diagram of the XX Z model in a 
longitudinal field [l9|. However, for non-zero anisotropy 
7, a gapless floating phase (FP) emerges between the SF 
and the AFM phases which is characterized by power- 
law decay of the correlations 



[131, 1201 121|. The transi- 
tion from the AFM to the FP is of the commensurate- 
incommensurate (C-IC) type whereas the transition be- 
tween the FP and SF phases is of the Bcrczinsky- 



Kosterlitz-Thoulcss (BKT) type. For A < -(1 + I7I) 
there is a first order transition aX h — Q between the 
two polarized phases. Finally, there is an Ising transition 
between the PP and the SF phases. 

Finite size 'phase diagram'. — The experimental real- 
ization of the Heisenberg model with cold atoms will in- 
evitably involve finite size effects due to the harmonic 
trapping potential. Within the local density approxima- 
tion, the trap renormalizes the couplings so that they 
become spatially dependent [22|, but this effect can be 
negligible if the orbitals are small compared to the length 
scale of the trap. In the regime of strong repulsion, the 
main effect of the trap is instead that it gives rise to 
"wedding cake" structures with Mott regions of integer 
filling. This effect was observed in the lowest band Bose- 
Hubbard model [l[ , and predicted theoretically to occur 
for antifcrromagnetic systems [23|. We therefore expect 
that for trapped systems, there will be finite Mott regions 
where the physics is controlled by the Heisenberg model 
as described in this paper. 

To examine finite size effects, we have performed an 
exact diagonalization of a chain with 18 spins with open 
boundary conditions. Figure [1] (b) displays the result- 
ing finite size 'phase diagram'. The colors correspond to 
different values of the total magnetization Af = J^ii^i) 
of the ground state. While the PP phase and the AMF 
phase are both clearly visible, the numerical results re- 
veal a step like structure of the magnetization in between 
the two phases. We attribute these steps in M to a 
devil's staircase structure of spin-density- waves (SDW). 
Consider a system of L sites, open boundary conditions 
and fixed values of 7 and A > 1 -I- |7|. Increasing the 
field from /i = 0, the staggered magnetization experi- 
ences modulations until the system enters the PP. The 
system passes through a succession of plateaus with dif- 
ferent values of total magnetization which build an in- 
complete devil's staircase l20[. As we see from Fig. [1] 
(b), the result is that it is only possible to give a numer- 
ical result for the PP-SF Ising transition. In particular, 
the C-IC and BKT transitions arc overshadowed by the 
transitions between SDW. Larger values of L allow for a 
larger number of plateaus and thus to a larger number of 
SDW with different periodicities. In the thermodynamic 
limit the staircase becomes complete and the changes in 
M become smooth. One then recovers the phase diagram 
of Fig. [T] (a). These transitions, between different SDW, 
are more pronounced for moderate systems sizes. For a 
typical experimental system with ~50 sites, for example, 
we estimate -^15 different SDW between the AFM and 
PP phases. 

Measurements and manipulations. — While time-of- 
flight measurements can reveal some of the phases [12l | , 
single-site addressing techniques 2J] will be much more 
powerful when extracting correlation functions. To ad- 
dress single orbital states or even perform spin rota- 
tions, one may borrow techniques developed for trapped 
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FIG. 2. (Color online) Different types of models are achieved 
by varying the relative tunneling strength and the relative 
orbital squeezing. The three different parameter regions are: 
(I) anti-ferromagnetic couplings in all spin components with 
A > J(l + I7I), (II) ferromagnetic or anti-ferromagnetic cou- 
plings in the ^-component and anti-ferromagnetic in the y- 
component with J(l + I7I) > |A|, and (III) same as in (II) 
but with |A| > J(H- I7I). The inset shows one example of the 
spin parameters J^x ~ (H- 7), Jyy — (1 — 7), and Jzz = A/ J 
for ty/tx = -0.1. 
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Making use of the p^ and Py orbitals symme- 
tries, stimulated Raman transitions can drive both side- 
band and carrier transitions for the chosen orbitals in 
the Lamb-Dicke regime. These transitions can be made 
so short that the system is essentially frozen during the 
operation. Driving side-band transitions in this way, spin 
rotations may be implemented. For example, a rotation 
around x is achieved by driving the red-sideband for both 
orbitals and in the large detuned case an adiabatic elimi- 
nation of the s-band results in the desired operation. Al- 
ternatively, Stark-shifting one of the p-orbitals results in 
a rotation around z. Since the spin operators do not com- 
mute, any rotation can be realized from these two opera- 
tions. Performing fluorescence on single orbital states by 
driving the carrier transition acts as measuring Sf . This 
combined with the above mentioned rotations makes it 
possible to measure the spin at any site in any direction. 
Tuning of couplings. — For a square optical lattice, 
we have Uxx ~ Uyy. Moreover, in the harmonic ap- 
proximation Uxy = Uxxl^i from which is follows that 
A < and 7 = —1/2. Thus, for a square lattice in 
the harmonic limit, interactions between the z— compo- 
nent of neighboring spins have ferromagnetic couplings, 
while the interactions between x— and between the y— 
components have anti-ferromagnetic ones. We now show 
how the relative strength and sign of the different cou- 
plings can be controlled by squeezing one of orbital 
states. Such squeezing can be accomplished by again 
driving the carrier transition of either of the two or- 
bitals, dispersively with a spatially dependent field. The 
shape of the drive can be chosen such that the result- 
ing Stark shift is weaker in the center of the sites, re- 



sulting in a narrowing of the orbital. To be specific, 
assume that the ratio of the harmonic length scales of 
the px and py orbitals in the y-direction, a^ is tuned. 
A straightforward calculation using harmonic oscillator 
functions yields a = Uxx /Uxy = 2^^/23(1 + o-^)^/^/cr and 
(3 = Uyy/Uxy = 2-3/23(1 + ^2)3/2^ rpj^g coupliug cou- 
stants now depend on a as A/J = 2<^(t^)~-'^/3/(a/3 — 
1) 4- 2t2'(P)"^"/(a/3 - 1) - [t^'/ty + ty/t'=)/2 and 7 = 
— 4/(q;/3 — 1). The inset in Fig. [5] displays the three cou- 
pling parameters as a function of a for |i^/i^| = 0.1. We 
sec that the relative size and even the sign of the cou- 
plings can be tuned by varying a. In particular, while 
Sy always has AFM couplings, they can be made both 
FM or AFM for 5":^ and S^. In the main part of Fig. [^ 
we sketch the different accessible models as a function 
of fy /t^ and a. This clearly demonstrates that one can 
realize a whole class of XYZ spin chains by using this 
method. 

Experimental realization. — In Ref. |ll| . the experi- 
mental realization of p-band bosons with a lifetime of 
several milliseconds was reported. With an average num- 
ber of approximately two atoms per site, the atoms could 
tunnel hundreds of times in the p-band before decay- 
ing. The main decay mechanism was atom collisions [9|. 
When there is only one atom per site, an increase in the 
lifetime is therefore expected and corresponds to a factor 
of five longer [ll| . Typical values of the couplings can be 
estimated from the overlap integrals of neighboring Wan- 
nier functions. Taking ®^Rb atoms and Xiat = 843n77i, 
we have J/i^j?, — 0.1 with Er = h k'^/2m, and and the 
characteristic tunneling time is t = h/JEn ~ 500 /iS. 
This is a few hundred of times smaller than the expected 
life-times [ll| which allows for the experimental explo- 
rations of our results. Furthermore, the spin correlations 
we discuss in this paper will emerge for temperatures 
ksT ^ J ^ i^/U. This is the usual temperature scale 
for spin models based on exchange couplings [7|. 

Conclusions. — We showed that the Mott regime of 
unit filling of bosonic atoms in the first excited band 
of an optical lattice realize the spin-1/2 XYZ quantum 
Heisenberg model. We then illustrated the rich physics 
of this model by examining the phase diagram of the 
ID case. Finite size effects relevant to the trapped sys- 
tems were shown to give rise to spin-density-waves and 
a devil's stair case structure. We proposed a method to 
control the strength and relative size of the spin couplings 
thereby demonstrating how one can realize a whole class 
of XYZ models. In this respect, the system is a promis- 
ing candidate for simulations of quantum magnetism. We 
finally discussed experimental issues related to the real- 
ization of this model. We end by noting that an exciting 
prospect is the 3D case, where we expect the three orbital 
model to result in even richer physics. 
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